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Abstract 



We present a new formalism which allows to derive the general Lagrangian dynamical 
equations for the motion of gravitating particles in a non-flat Friedmann universe with 
arbitrary density parameter fl and no cosmological constant. We treat the set of particles 
as a Newtonian collisionless fluid. The non-linear dynamical evolution of the fluid element 
trajectories is then described up to the third-order expansion in Lagrangian coordinates. 
This work generalizes recent investigations carried out by Bouchet et al. (1992) and Buchert 
(1994). 

Subject headings: Galaxies: clustering - large-scale structure of the Universe 

1 Introduction 

The most widely accepted hypothesis about the formation of large-scale structures is that 
galaxies and clusters of galaxies formed by gravitational collapse around primordial slightly 
overdense mass fluctuations in the Universe. One way to link these initial conditions to the 
final mass distribution is to attempt a solution of the equations for the matter fluctuation 
field 5 and the peculiar velocity field v, namely the Euler equation, the continuity equation 
and the Poisson equation. 

Linear theory easily provides solutions even for cosmological models with arbitrary 
density parameter Q (e.g. Peebles 1980), but when the perturbation amplitudes approach 
unity, non-linear gravitational effects become very important. It is undoubtedly impossible 
to follow the non-linear stage of the gravitational collapse in an exact analytical way, 
and one is forced to use approximation techniques such as perturbation theory where, for 
instance, the systematic expansion of the density solution is obtained by writing 5 = J2 
where 5^ = 0(5^ n ), 5^ corresponding to the linear solution (Fry 1984; Goroff et al. 1986; 
Juszkiewicz, Bouchet & Colombi 1993; Catelan & Moscardini 1994a,b; Bernardeau 1992, 
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1994; Catelan et al. 1995). 

Usually, these analyses are intrinsically Eulerian, the fundamental quantities being the 
density and velocity fields evaluated at the (comoving) Eulerian coordinate x; also, with the 
exception of Bernardeau (1992, 1994) and Catelan et al. (1995), they have been confined 
within the limits of the Einstein-de Sitter cosmology, mainly because the condition Q = 1 
enormously eases the investigations of solutions of the dynamical equations for 5 and v: 
apart from theoretical expectations, there is however no definitive observational evidence 
that our Universe is really flat (see e.g. Peebles 1991; Coles & Ellis 1994). 

It has been recognized that the problem of giving an analytical description of the non- 
linear process of gravitational clustering simplifies when formulated in terms of Lagrangian 
coordinates rather than the standard Eulerian ones: Zel'dovich (1970a, b) first proposed to 
approximately describe the weakly non-linear regime of density evolution in terms of the 
departure from to the Lagrangian (initial) positions of the fluid elements. The Zel'dovich 
approximation is now widely used in cosmology, showing also to be extremely useful in 
reconstruction methods of initial conditions from velocity data (e.g. Nusser & Dekel 1992). 
However, only recently it has been fully understood that the whole dynamics of gravitational 
clustering may be suitably described in terms of the displacement field S, which turns out, in 
the Lagrangian approach, to be the only underlying fundamental field. Buchert (1989, 1992) 
indeed derives the exact equations governing the evolution of the displacement S (therefore 
of the density and velocity fields): since, as in the Eulerian case, it is impossible to work out 
the general solution S, a perturbative approach is again introduced. The key novelty with 
respect to the Eulerian approach is that one searches for solutions of perturbed trajectories 
about the linear (initial) displacement : S = J2 where S^ n ^ = o(S^) (see Moutarde 
et al. 1991). The important point is that a slight perturbation of the Lagrangian particle 
paths carries a large amount of non-linear information about the corresponding Eulerian 
evolved observables, since the Lagrangian picture is intrinsically non-linear in the density 
field. 
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Solutions up to the third-order Lagrangian approximation have been obtained (the 
first-order solution corresponding to the Zel'dovich approximation), although limited to 
the case of an Einstein-de Sitter model (Buchert & Ehlers 1993; Buchert 1994). The 
higher accuracy of Lagrangian perturbative methods as compared to other currently studied 
approximation ansatzs [such as the frozen flow (Matarrese et al. 1992) and the linear 
potential approximations (Brainerd, Scherrer & Villumsen 1993; Bagla & Padmanabhan 
1994)] is discussed in Munshi & Starobinsky (1994), Bernardeau et al. (1994) and Munshi, 
Sahni & Starobinsky (1994), again in the framework of an Einstein-de Sitter cosmology. 
Comparisons with iV-body simulations in the fully developed non-linear clustering are 
displayed in Moutarde et al. (1991), Coles, Melott & Shandarin (1993) and Melott et 
al. (1994). 

The second-order Lagrangian solution for generic non-flat Friedmann models has been 
derived by Bouchet et al. (1992; hereafter BJCP), where particular emphasis on the con- 
nection with the Eulerian formulation has been given. Further attempts to extend the 
Lagrangian formalism to models with arbitrary density parameter may be found in Gra- 
mann (1993) and Lachieze-Rey (1993b), which, however, lead to not completely correct 
conclusions. Matarrese, Pantano & Saez (1994a,b), developed a relativistic Lagrangian 
treatment of the non-linear dynamics of an irrotational collisionless fluid, which reduced 
to the standard Newtonian approach on sub-horizon scales, but is also suitable for the 
description of perturbations on super-horizon scales. 

In this work, we present a new Lagrangian formalism which enables one to easily derive 
the exact dynamical equations governing a pressureless Newtonian gravitating fluid in an 
expanding universe with arbitrary density parameter Q and no cosmo logical constant. We 
then describe the non-linear evolution of perturbations up to the third-order Lagrangian 
approximation. The first-order solution corresponds, of course, to the Zel'dovich approxi- 
mation, and the second-order one to the BJCP approximation. The third-order solution 
is then derived in detail for arbitrary initial conditions. It consists of three independent 
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(growing) modes, two being purely longitudinal and one purely transversal, in such a way 
that a fluid which is irrotational in Eulerian space is not so in the corresponding Lagrangian 
space. This problem was first addressed in Buchert (1994), where the consequences of re- 
quiring irrotationality of the fluid motion in Lagrangian space on the initial conditions are 
explored in detail. 

The layout of this paper is as follows. In section 2 the Lagrangian approach is reviewed. 
In section 3, the general Lagrangian equations describing the dynamical evolution of a col- 
lisionless Newtonian self-gravitating fluid in an expanding universe with arbitrary density 
parameter Q are obtained. In section 4, after rederiving in our formalism the Zel'dovich 
and the BJCP approximations, we work out the third-order Lagrangian approximation. 
Our conclusions are presented in section 5. To avoid an excess of mathematical contortions 
in the text, four technical appendices are given. 

2 Lagrangian Formulation 

Let us consider a Newtonian pressureless self-gravitating fluid embedded in an expanding 
universe with arbitrary density parameter fl. We assume that the cosmological constant is 
exactly zero. We have in mind that such a fluid mimics the behavior of matter on scales 
smaller than the horizon scale and that, around its primordial density perturbations 5, the 
present luminous objects like galaxies or clusters of galaxies started to grow according to 
a gravitational instability process. We indicate by x the comoving Eulerian coordinates, 
from which physical distances may be obtained according to the law r = a(t)x, a(t) being 
the expansion scale factor and t the standard cosmic time. 

According to the Lagrangian point of view, the path of each fluid element is followed 
during its evolution. A proper observer may label each neighbouring fluid particle by e.g. 
its initial (comoving) coordinate, say q = x Q . At time t, the (Eulerian) position of the 
q-particle will be 

x = x(q,*)=x L . (1) 
5 



Here the only independent variables are the labels q (apart from the time t) which therefore 
play the role of spatial coordinates in the Lagrangian {q}-space. The vector x, which in the 
Eulerian picture is an independent variable, is now introduced as a new real dynamical field: 
in the Lagrangian {q}-space, the trajectories of the mass elements are fully described by the 
dynamical maps x(q, t), starting from the initial positions q. The definition (1) implicitly 
assumes that there is a one-to-one correspondence between the Eulerian coordinate x and 
the Lagrangian coordinate q: this is certainly the case for a cold non collisional fluid, at least 
until the stage of caustic formation (see e.g. Shandarin & Zel'dovich 1989). Mathematically, 
this is equivalent to the statement that the determinant J of the Jacobian of the map 
q — > x(q, t) is non-singular, 



the map x(q, t) being thus reversible to q(x, t). Obviously, during the highly non-linear 
evolution, many particles coming from very different original positions will tend to arrive 
at the same Eulerian place: infinite-density regions (caustics) will form in Eulerian space 
and the map from Lagrangian to Eulerian space becomes singular (Shandarin & Zel'dovich 
1989; see also the discussion in Kofman et al. 1994). This caustic formation process enor- 
mously limits the predictive power of any (perturbative) Eulerian method, in that the 
main requirement for it to work is the very restrictive condition |<5| <C 1. This is possibly 
the strongest reason for preferring the alternative Lagrangian picture: in this description 
indeed the density is not a dynamical variable and is fully integrated. 

Following Zel'dovich (1970a, b), all Eulerian fields can be represented in terms of the 
only dynamical field x^ - or its derivatives. Thus, the peculiar velocity u, the peculiar 
acceleration f and the density p are respectively given by 




(2) 



<fo(q, t) 
dt 



u(q, t)a{t) 



-i 



(3) 



x(q,t) 
dt 2 



[f(q,t)-2ifu(q,f)]a(t) 



i 



(4) 



l + 5[x(q,t),t] = J(q,t) 



-l 



(5) 
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the quantity H being the Hubble constant and 5 = (p — pb)/pb the density fluctuation about 
the density background value Pbif). The operator d/dt is the usual convective Lagrangian 
time derivative which follows the mass element, ^ = J^| q = ^|x + ^- V x • Recall that the 
operator ^ does not commute with the Eulerian nabla operator V x . 

Equation (4) corresponds to the (comoving) Euler equation and eq. (5) to the continuity 
equation for the system. Furthermore, the equations (3) and (4) - governing the evolution of 
the field = x(q, t) - may be thought as definitions of the velocity u and the acceleration 
f . The particular form of the mass conservation (4) derives from the fact that the observer 
follows the mass elements during their motion and that the initial matter distribution is 
assumed homogeneous. 

The gravitational interaction among the mass elements is then introduced, requiring 
that the peculiar acceleration field f is induced by the density fluctuations 5 through the 
Poisson equation, V • f = —4nGapb 5. For our purposes, however, the best form of the 
Poisson equation is the one obtained introducing the time variable (Doroshkevich, Ryabenki 
& Shandarin 1973; Shandarin 1980) 

r = V^k(l-Q)-^ 2 , (6) 

where k — —1 for open universes and k = 1 for closed universes. The case Q = 1 (k — 0) is 
a singular point for the transformation (6), and, in that case we take r = t -1 / 3 . As we shall 
see below, the main advantage in using r instead of the standard cosmic time t is that the 
dynamical equations for the higher-order Lagrangian modes are considerably simplified, 
since the rescaled peculiar acceleration g is now simply given by 

Finally, gravity is introduced through the equation 

V x -g = -a(r)(J, (8) 

where the quantity a(r) = 6/(r 2 + k) (Shandarin 1980). The time variable r has been 
recently used also in BJCP and Lachieze-Rey (1993a,b). 



Alternatively, if the original f is expressed in terms of the gravitational potential $, i.e. 
f = — a _1 V$, then, because of (8), 

V^ = a(r)5, (9) 

once the rescaled potential = (2a/VL a a 2 ) 2 <3> such that g = — V<p has been introduced. 

Finally, if one assumes irrotational motions in Eulerian space, which is a plausible hy- 
pothesis for a collisionless (cold) dust, then the "irrotationality" condition may be expressed 
according to the relation 

V x Au = 0. (10) 

Summarizing: the Lagrangian system of equations governing a gravitating collisionless 
fluid is given by eqs.(3) and (4) - the Euler equation - for the field x(q, t), the mass 
conservation relation (5) and the Poisson equation. 



3 Lagrangian Equations for the Displacement 

To describe the departure of the mass elements from the initial positions q one usually 
introduces the displacement vector S such that 

x(q,r) = q+S(q,r). (11) 

It is clear from this definition that the motion of the fluid elements may be completely 
described in terms of the displacement S, since the latter fully characterizes the map (1) 
between the Eulerian and the Lagrangian coordinates. In terms of S, the Euler equation 
and the continuity equation may be written, respectively, as 

d 2 S(q, t) r / x i ,,~n 

™ ' = g[x(q, r),r] , (12) 

l + 5[x(q,r),r] =det(/ + 5)- 1 . (13) 

Here I = diag(l, 1, 1) and S is a 3 x 3 matrix whose elements are S a/ 3 = dS a /dq/3, also 
called the deformation tensor. In general the deformation tensor is not symmetric, i.e. 



S a /3 7^ S/3 a : S a p is symmetric iff the displacement S is potential in the Lagrangian space. 
The Eulerian "irrotationality" condition (10) - if taken into account - and the Poisson 
equation (8) may be written as 

V x AS = 0, (14) 

V x -S = -a(r)<J[x(q,r),r] , (15) 

the velocity field being defined by the relation dx/dr = a(r) u(x, r). The second equation 
clearly shows how the trajectories S are deformed during the time evolution by the density 
perturbation <5[x(q, r), r]. The dot indicates the operator d/dr. However, we cannot say 
that the previous equations are completely written in the Lagrangian {q}- space: indeed 
the operator V x does not act on the Lagrangian coordinates q. The exact way in which 
the differentiation with respect to the Eulerian position x is translated - through the map 
(11) - into differentiation with respect to the Lagrangian position q is displayed in the 
following relation: 



J 



d 



d 



dxp L ;i + v.s)^-^ + ^j-, (i6) 

where now V = V q . The symbol 5 a p indicates the Kronecker tensor and the quantity S%p 
is an element of the cofactor matrix S . Summation over repeated Greek indices (where 
a = 1, 2, 3) is understood. We derive the above relation in detail in Appendix A. 
The Newtonian Lagrangian equations for the collisionless fluid finally become 



and 



1 + V • S) V " Sfr + Sg, S Ja = , (17) 



1 + V • S) 5 al3 - S af3 + S^\ S Pa = a(r)[J(q, r) - 1] , (18) 

where e a/ g 7 is the totally antisymmetric Levi-Civita tensor of rank three, e i2 3 = 1. 

We call the latter equation the Lagrangian Poisson equation. The irrotationality con- 
dition in Lagrangian space (17) and the Lagrangian Poisson equation (18) dynamically 
constrain (up to a constant vector) the field S. They are the (closed set of) general dy- 
namical equations for the displacement vector S describing the motion of a collisionless 
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fluid in the Lagrangian {q}- space, embedded in an arbitrary non-flat universe and subject 
to the Newtonian gravitational influence of the mass fluctuations J -1 — 1. We stress that 
eqs.(17) and (18) are non-linear and non-local in the displacement S (see the discussion in 
Kofman & Pogosyan 1994). Analogous equations, although using a very different tensorial 
notation, are analyzed in Buchert (1989). 

Some remarks are appropriate. Eq.(17) does not mean that the motion is potential in 
the Lagrangian space, since this would correspond to the condition e a py Spy = 0. On the 
contrary, we can surely state from (17) that the motion is vortical in Lagrangian space and 
that, however, the departure from the irrotationality of the Lagrangian peculiar velocity 
u is gravitationally induced only at higher order in the displacement S, as may be clearly 
seen if one writes e^S^ = — e Q( g 7 ^^[(V • S) 5 la — S 1C + S^ a ]. This implies that, at 
least to the first order in S, the condition for irrotationality in Eulerian space means that 
the Lagrangian motion is also potential, as first noted by Zel'dovich (1970a, b). The 
irrotationality problem in Lagrangian coordinates has been recently addressed by Buchert 
(1992, 1994). We will discuss again it below for a general nonflat model. Considering now 
the Lagrangian Poisson equation, we stress the fact that only the left hand side of eq.(18) 
is the 'dynamical' part, containing the time derivatives of the field S: note that a term 
proportional to S is absent; furthermore, the Lagrangian Poisson equation is intrinsically 
non-linear in the density field, unlike the Eulerian Poisson equation, as it may be seen from 
the relation J — 1 = — 1, where 5 is fully integrated. 

It is remarkable the fact that the equations (17) and (18) hold for a generic non-flat 
model, in that Q - although fundamental - is a very poorly known cosmological parameter 
and, apart from theoretical expectations (essentially due to the implications of inflation), 
there is no definitive observational evidence that our Universe is flat (Peebles 1991; Coles 
& Ellis 1994). Furthermore, eqs.(17) and (18) are manifestly comoving, which is also 
convenient, because the overall expansion is usually subtracted in analytical or numerical 
studies of departures from the mean homogeneous and isotropic universe. In Appendix A, 
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we give more compact expressions for eqs.(17) and (18). 

The irrotationality condition and the Lagrangian Poisson equation are exact equations 
in the Lagrangian description. It is undoubtedly very difficult to solve them in a rigorous 
way. A possible alternative is to seek for approximate solutions: the standard technique 
is to expand the trajectory S in a perturbative series, the leading term being the linear 
displacement which corresponds indeed to the Zel'dovich approximation (see Moutarde et 
al. 1991). To approximate the Lagrangian Poisson equation implies that the gravitational 
interaction among the particles of the fluid is described only approximately. 

4 Lagrangian Perturbative Approximation: Higher- 
Order Solutions 

We now solve the dynamical equations for the displacements S according to the following 
Lagrangian perturbative prescription: 

S(q, r) = D{r) S«(q) + E(r) S< 2 >(q) + F(r) S< 3 >(q) + • • • . (19) 

Here S^^(q) corresponds to the first-order approximation, S^(q) to the second-order 
approximation, and so on: the dynamics of the evolution constrains in general both the 
temporal dependence as described by the functions D, E, F, and the spatial displace- 
ments S (n) (q). 

Recently Gramann (1993), to calculate an analytical expression relating the density 
to the velocity in a self-gravitating system, used a (second-order) perturbative expansion 
similar to the previous one, but with the assumption that E = D 2 : we will discuss below 
why this is a very restrictive hypothesis. Furthermore, we anticipate here - and we will 
demonstrate in Section 4.3 - that the third-order solution actually corresponds to 
three modes: we maintain here the expression as given in eq.(19) for the sake of simplicity. 

We emphasize that the first- and the second-order perturbative solutions (and, as it will 
be shown, each mode of the third-order solution) are explicitly written in (19) as separable 
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with respect to the temporal and spatial coordinates. This is not an assumption, being just 
a property of the perturbative Lagrangian description. Indeed, this can be demonstrated 
to be just a direct consequence of the dynamical Lagrangian equations: while this turns 
out to be trivial for the first-order solution, it is not so e.g. for the second-order one. 
In Appendix B, it is shown that the separable solution E(r)S^ (q) is the most general 
second-order Lagrangian perturbative solution. 

Conversely, in Catelan et al. (1995) the non-separability of the corresponding higher- 
order Eulerian perturbative solutions is thoroughly analyzed: the Eulerian perturbative 
solutions factorize in space and time only in the very special case of the flat universe, and 
not in a generic Friedmann model; intriguingly, it is also shown that the non-separability of 
the Eulerian solutions is fully consistent with the separability of the Lagrangian solutions, 
at least explicitly up to the second-order. 

Because of their convenience in deriving the subsequent equations, we define the follow- 
ing scalars: 

Mi (A) = V ■ A = A aa , (20) 

/i 2 (A, B) = X - [A aa Bf3p - A a(3 B Pa ] , (21) 

/! 2 (A) = /i 2 (A,A), (22) 

/i 3 (A) = det(A al3 ) , (23) 

where A and B are generic functions of the Lagrangian coordinate q and spatial derivatives 
are with respect to the variable q. Note that the functions fj,i, yU 2 and /i 3 are linear, quadratic 
and cubic in their arguments, respectively. Furthermore, the following expression for the 
Jacobian determinant J holds: 

J(q, r) = 1 + MS) + // 2 (S) + // 3 (S) . (24) 

This is an exact relation for the displacement and it is cubic in S, justifying investigation 
of the third-order solution. According to eqs.(12) and (19), the Lagrangian acceleration 
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may be perturbatively expressed as 

g[x(q, r),r]= D(r) S«(q) + E{r) S< 2 >(q) + F(r) S< 3 >(q) + • • • , 



(25) 



explicitly up to the third-order term: this expression corresponds to the Euler equation. 
We now solve the Lagrangian Poisson equation order-by-order. 

4.1 First— Order Solution: Zel'dovich Approximation 

We can easily find the first-order approximation truncating eq.(18) accordingly to the linear 
terms, S aa = «(t) /ii(S) , to obtain, in terms of the displacement 

D(r) = a(r) D(r) ^) . (26) 

Given the definition of /ii, we immediately get 

D(t) - a(r) D{t) = . (27) 

The two linearly independent solutions coincide with the growing and decreasing modes of 
the linear density field: hereafter we will consider only the growing mode D + = D, since any 
perturbative approach is consistently applicable to the (mildly) non-linear regime, when 
the decreasing modes have already become negligible. We report here the linear growing 
solution 

D(r) = l + 3(r 2 -l] 



1 + rim 



It -I 



r + 1 



(28) 



where we remind that, in the case k = 1, ln(^) = 2iarctang(ir). Eq.(27) thus corre- 
sponds to the Zel'dovich approximation. Note that the linear regime does not constrain 
the vector at all, whose particular form has to be ascribed to the chosen initial con- 
ditions. Furthermore, the displacement is potential in Lagrangian space, since the 
irrotationality condition (17) in the linear regime reduces to 

e a ^S^ = 0, (29) 
13 



as already known (Zel'dovich 1970a, b). One can thus define a potential tp^(q) such 
that S^(q) = V^^^(q), with ip^ the velocity potential in the Zel'dovich approximation 
(see the discussion in Kofman 1991). As a consequence, the linear deformation tensor is 
symmetric. The particular growing solution (28) reduces to D — t~ 2 in the case of an 
Einstein-de Sitter model, since a = 6t~ 2 . 

4.2 Second— Order Solution: BJCP approximation 

Retaining only the quadratic terms in the Lagrangian Poisson equation and introducing 
the ansatz S = DS^ + ES^ 2 \ one gets to second-order 

EMS^ + DZM^) 2 -^^ , 

from which, using the first-order results, 

E(t) - a(r)E(r)] ^(8®) = -a(r) D{r f /i 2 (S«) . (30) 

This is a separable differential equation leading to the system of equations (cf. BJCP 1992) 



E{t) - a{r)E{r) = -a(r) D{r) 2 , 
/ii(S( 2 ')=/i 2 (SW). 



(31) 



We stress that, unlike the first-order case, the second-order approximation constrains both 
the time and the spatial dependence of the solution. The solution E(r) of the temporal 
equation has been found by BJCP (1992). Its generality is discussed in Appendix B. We 
report here its explicit expression: 

2' 



E(r) = - 1 -- 9 -(r 2 -l){ 



r + (r 2 - l)lm 



T+ 1 



(32) 



In BJCP, this solution has been applied to describe the dependence of the skewness of the 
unfiltered density field on the density parameter f2, in that, near Q, — 1, E « — |f2~ 2 / 63 .D 2 . 
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For us it will be useful, because we will be able to express the third-order solutions in 
terms of the lower-order results. 

The solution of the second-order spatial equation in (31) may be written as 

S ( 2 ) = - [s« (V • S«) - (S« • V) S«] + R( 2 ) , (33) 

where R^(q) is any divergence-free vector such that V A S*^ = 0: if the first-order 
Lagrangian motion is assumed potential, this is indeed the result one obtains from the 
irrotationality equation (17) once only the second-order terms are retained 

E{r) (V A S^) a = -D{r)D{r) e a „ S$S$ = , (34) 

where the last equality follows from the fact that the tensor sfjs^j is symmetric. Thus, we 
have to conclude that the gravitational evolution does not induce vorticity in the second- 
order Lagrangian motion at all, whatever the initial conditions are: note that VAS^ = is 
a purely spatial relation. It appears clear from (34) that the displacement is potential 
once S^ 1 ) is assumed to be potential too, and thus the second-order deformation tensor is 
symmetric: a potential may be introduced, from which one obtains S^(q) = V^ 2 )(q). 
An useful expression of the Fourier transform of ip( 2 \ W 2 \ in terms of the first-order 
potential may be found in Appendix C: it will clearly result that an explicit expression 
of the gauge-dependent vector R^ 2 -* is, in practice, unnecessary: we nevertheless emphasize 
that the expression in (33) is the only compatible with the irrotationality condition (34). 
A restricted class of second-order irrotational solutions, for which additional constraints 
on the initial conditions have to be fulfilled, is discussed in Buchert & Ehlers (1993) and 
Buchert (1994); the same constraints also allow the construction of local forms, as debated 
in Buchert (1994), again in the context of the Einstein-de Sitter cosmology. The solution 
in (32) reduces to the simple form E = — | t~ 4 oc D 2 in the case of the flat model. 
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4.3 Third-Order Solution 



Inserting the ansatz (19) in the Lagrangian Poisson equation, after some algebra one obtains 
the third-order expression 



F/ii(S (3) ) + 2(DE + DE) /i 2 (S (1) , S (2) ) + 3DD 2 ^(S 



a(r) If^S^) + 2DEfi 2 (S^, S< 2 >) + D 3 /x 3 (S«) 



where we used the relation 



>2 o(l)C c (l) 



On the light of the lower-order results, the former expression becomes 



(35) 



F(r) - a(r)F(r)] ^(S®) 

= -2a(r)£>(r) [^(r) - D{rf] ^(S«, S< 2 >) - 2a(r) J D(r)-V 3 (S (1) ) . (36) 

Note that the mixed invariant /i2 couples the first- and second-order solutions and S^ 2 \ 
while the last term // 3 is cubic in the argument This fact forces us, as suggested by 
Buchert (1994), to split the third-order displacement into two parts, one resulting from 
the interaction among the linear perturbations, the second from the interaction between 
the first- and second-order perturbations: 



S (3) (q) = Sf(q) + Sr ) (q). 
According to this ansatz, the dynamical equation (36) splits into the systems 

F a (r)-a(r)F a (r) = -2a(r)D(r) 3 , 
^(Si 3) )=ft(S«), 



(37) 



(38) 



and 
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F b {r) - a{r)F b {r) = -2a(r) D(r) [E{r) - D{rf \ , 

(39) 

Ms[ 3) ) = Ms« s( 2 )). 

Note that, since now may be considered a potential field, /i 2 (S^\ S^) = /i 2 (S^ 2 \ S^). 

It clearly appears that the assumption E = D 2 [see Gramann (1993)] implies a neglect 
of the interaction mode between the first- and second-order perturbations, in that E = 
D 2 ==>- F b = 0. This is in general a very restrictive hypothesis and possibly only (high- 
resolution) iV-body simulations might quantify the real loss of accuracy at the third-order 
level. 

Furthermore, Lachieze-Rey (1993b) argues that one formal third-order solution exists 
corresponding to the assumption that the deformation tensor <S(q, r) remains proportional, 
at each Lagrangian position, to its initial value, being only multiplied by a scalar Lagrangian 
growth factor. This is surely true, but the problem is now that such a hypothesis notably 
restricts the form of the spatial part of the third-order solution: it is not difficult indeed 
to show that the Lachieze-Rey solution corresponds to assuming that, in our notation, 
^ 3 (SW) = /i 2 (S (1) ,S (2) ) and, after that, ^i(S (3) ) = yU 3 (S (1) ). Of course these conditions 
depend also on the particular chosen initial configurations, but it is not clear to which type 
of realistic physical situation they may be applied. 

After these comments, we give now the third-order solutions by quadrature: 

F a (r) = -2D(r) f dr x D^)' 2 f' dr 2 a(r 2 ) D(r 2 ) 4 , (40) 

and 

F b (r) = -2 D{t) f dn Din)- 2 P dr 2 a(r 2 ) D(r 2 ) 2 [e(t 2 ) - D(r 2 ) 2 } . (41) 

In this way the two solutions F a {r) and F h {r) are given in terms of the lower-order solutions 
D{t) and E{r): explicit versions of (40) and (41) are cumbersome, while approximate 
expressions may be recovered in the limit e.g. Q — > 1 (Catelan 1995). In the case of an 
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Einstein-de Sitter model, we find F a = — |r 6 and F b = 6 oc D 3 : the results in Buchert 
(1994) are particular cases of the general solutions (40) and (41). 

A solution S< 3 > of (38) may be obtained recalling that, from (35), S^f Sfl = 3// 3 (S«); 
thus, since d a S^ C = 0, one gets d a (S^ — \S^ C S^^j = and finally, for each component, 

S^ = ls^Sjp+R^l. (42) 

Here R® (q) is again a divergence-free vector such that V A = 0. To understand why 
is a potential vector it is indeed enough to isolate, in the Lagrangian irrotationality 
condition (17), all the possible vortical terms induced by the interactions among the first- 
order displacements: it is not difficult to see that the only term 

e a ^sSs$ c = 0, (43) 

where the equality is justified by the assumption that the linear displacement is potential. 
Another way to describe the irrotationality of is to note the that third-order interac- 
tions like Si]) SiV SiV can be only symmetric in SiV and not antisymmetric (here * indicate 
a generic Greek index), as any possible vortical-like interaction would be. 

In a similar fashion, one can derive a solution S b from (39), assuming in particular that 

(3) 

S b is induced by symmetric couplings between the first- and the second-order solutions. 
It can be written as 



si 3) = 



1 r 



S«(V • S< 2 >) - (S« • V)S< 2 > + S< 2 >(V • S«) - (S< 2 > • V)S«1 + R[ 3) , (44) 



where again H b (q) is a divergence-free vector such that V A S b = 0: in fact, symmetric 
interactions between the first- and second-order modes cannot be vortical. 

The solutions and Sf \ being longitudinal, may be written in terms of the respective 
potentials ip^ and ip b 3 \ namely S^(q) = V^ 3 )(q) and S^(q) = VV^(q). In Appendix 
C, the expressions of the Fourier components ip^ and ip b 3 ^ as functions of the Zel'dovich 
potential ■ip^ are given. The significance of the vectors R^ and R^ is perfectly equivalent 
to that of the vector R( 2 ) for the second-order solution S^ 2 ^. 
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One now has to make the following inquiry: does any antisymmetric coupling between 
the Zel'dovich and the BJCP solutions exist at the third-order level? If it exists, then 
a third-order transverse component T arises in the Lagrangian motion of the fluid mass 
elements, and this is the only vortical component which can be induced in the framework 
of the third-order approximation. Note that such a type of component cannot modify the 
amplitude of the density fluctuation field. 

Let us suppose therefore that such a transverse mode actually exists: if the original 
ansatz (19) is now improved to 

S(q, r) = S,| + S ± = S,| + F c (r) T(q) , (45) 

where we recall that Sy = £>S (1) + £S (2) + F a S^ + F b sf \ it follows that, in terms of T, 
the Lagrangian irrotationality condition becomes 



F c (r)(VAT) a = D(t)E(t) - D(t)E(t) 



ta/3 7 iJp c J al ■ 



(46) 



We derive eq.(46) explicitly in Appendix D. Again, this equation may be split into a tem- 
poral and a spatial part: 



F c (t) = -o(r) D(rf , 

(VAT) a = ^s£s(;. 



(47) 



The growth factor F c (t) may be explicitly written by quadrature 

F c (t) = -j T dn f 1 dr 2 a(r 2 ) D{r 2 f . (48) 

This solution reduces to the simple F c = —^t~ 6 in the case of an Einstein-de Sitter universe. 

From the second equation in (47), we note - consistently - that the antisymmetric part 
of the third-order interaction between the modes and S^ 2 ^ generates the transverse 
component 



T(q) 



;S (1) • V)S (2) - (S (2) • V)S (1) 1 + Vv?(q) . 



(49) 
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In Appendix C the Fourier components T a are explicitly given. The term V<y2 is such that 
V • T = or, equivalently, the transverse component T can now be described as the curl of 
a vector potential, T = V A A, which cannot in general be removed by a suitable gauge- 
fixing of the initial conditions: one can say that the Lagrangian motion is no longer purely 
potential from third-order onward. 

This fact was first discovered by Buchert (1994) within the Einstein-de Sitter model, 
and therefore seems a universal feature of the Lagrangian motion, independently of the 
value of the density parameter Q (i.e. of the underlying model of universe). However, 
it is important to stress again that the displacements and S[ 3 -* are in general purely 
potential just because they originate through symmetric couplings among the lower-order 
longitudinal perturbations; also, the transverse displacement T is purely vortical because 
it originates through the antisymmetric coupling between the first- and second-order so- 
lutions, independently of the peculiar initial conditions one is picking out. We summarize 
the results of this work in the next section. 

5 Summary and Conclusions 

In this work we studied, in the Lagrangian description, the behaviour of a self-gravitating 
collisionless fluid embedded in a generic non-flat expanding universe. We assumed the 
Newtonian limit and zero cosmological constant. We have mainly focussed on the formal 
aspects of this kind of analysis, also proposing a new Lagrangian formalism. In particu- 
lar, we derived the general equations, in the forms given in (17) and (18), governing the 
Lagrangian motion of the mass elements in a universe with arbitrary density parameter 
Q and directly in comoving coordinates. We consider this kind of derivation one of the 
main results of this work. An alternative derivation and notation may be found in Buchert 
(1989). Since it is in practice very difficult to find exact solutions of the foregoing general 
equations, we solved them according to a Lagrangian perturbative approach, namely we 
sought approximate solutions S = J2 about the linear displacement as pioneered 
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by Moutarde et al. (1991). Our formalism enables one to easily recover the already known 
lower-order solutions, in particular the linear Zel'dovich approximation (Zel'dovich 1970a, 
b) and the second-order BJCP approximation (Bouchet et al. 1992). Then, we explicitly 
worked out the third-order solution, generalizing to an arbitrary Friedmann universe the 
recent results obtained by Buchert (1994) in the context of the Einstein-de Sitter cosmol- 
ogy. The question of the irrotationality in the Lagrangian space has been analyzed too. 
In particular, we found that the spatial solutions S^, are purely potential for 

any acceptable initial conditions - the only underlying hypothesis being that the linear 
displacement is longitudinal - as one finds if one carefully applies, order-by-order, the 
irrotationality condition (17). Similarly, the third-order transverse component T is purely 
vortical for any realistic initial conditions. Obviously, the reason of the existence of the 
transverse component T is that the transformation from Eulerian to Lagrangian coordi- 
nates x — > q is in general a non-Galilean transformation, and fictitious forces are induced 
in the Lagrangian description: no new physics can indeed appear. 

To avoid an exceeding proliferation of formulae in the main text of the article, we 
performed in Appendix C the Fourier analysis of the higher-order Lagrangian solutions: 
such a results can be important for numerical and practical applications. Specifically, 
we found that the third-order dynamics is fully described in terms of the tetra-potential 
^ = A 1 , A 2 ), since only two of the three components of the vector potential A 

are physically significative. 

Our results offer a tool with which to follow the dynamics of the formation of the struc- 
tures in the universe, as resulting from the non-linear gravitational instability. 
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Note added in proof After the submission of this article, Thomas Buchert addressed 
our attention on the Appendix A in Buchert (1989) where, employing the same temporal 
parameter in (6), an alternative expression of the equation (18) is given. The temporal 
coordinate (6) has been recently used also in Bharadwaj, S. 1994, ApJ, 428, 419 and 
Bouchet, F.R., et al. 1994, A&A, submitted. 
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Appendix A 

In the first appendix we explicitly show how to obtain the general Lagrangian equations 
(17) and (18). Let us start by deriving the Lagrangian Poisson equation: from eq.(15), this 
can be written in the form 

dS a (q,r) 



8x r 



a(r)[J(q,r)-l], 



(50) 



where J is the determinant of the Jacobian (§^), and x = q + S(q, r). To derive (18), all 
we need is to show that the following relation holds 



J 



d 

dx a 



1 + S-y-y) 5[3a — S[3a + S 



d 



This indeed corresponds to the implicit expression 

d dqp d ( dx a \ 1 d 



dx a dx a dq p \ dqp J dqp ' 
The quantities (§^) are elements of the matrix = (f^)- Since, explicitly, 



(51) 



(52) 



dq, 



1 + Sn 


Sn 


5*13 


5*21 


1 + S 22 


5*23 


S31 


S32 


i + s. 



33 



(53) 



we obtain, in matrix notation, 

(|)" = ){(l + V-S)/-5 + ^}, (54) 

where I = diag(l, 1, 1) is the identity, S = (Sap) is the deformation matrix and S c is the 
cofactor matrix we usefully write below: 



S 



c 



S22S33 — S32S23 



5*31 5*23 — S21S33 



S21S32 — S31S22 



S32S13 — 5'l25'33 



S11S33 — 5*31 •S'is 



5*31 Sn — S11S32 



S12S23 — S22S13 



S21S13 — SIIS23 



S11S22 — S21S12 



(55) 
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With these results one now immediately recovers the Lagrangian Poisson equation from 
the original eq.(15). 

In exactly the same way one can obtain the irrotationality condition in the Lagrangian 
space from the Eulerian eq.(10): 

dxp 

Recalling eqs.(ll), (12) and (52), this becomes 



dx,yd Sl _ o (5?) 



, dq a J dq a 

If now the relation (51) is inserted in the last equation, we get the general equation (17). 

To conclude we propose here even more compact forms for the irrotationality condition 
(17) and the Lagrangian Poisson equation (18): in fact, defining the cofactor element of 

*% = Jx$=[(l + V- S) 5 aP - S aP + SgJ , (58) 
where x~p = (§^) , we m ay finally write 

e Q( 3 7 Xp a Xy a = , (59) 

x%x Pa = a{r)[J - 1] , (60) 

respectively. We have nevertheless preferred to retain less compact and elegant versions in 
the main text, because the dependence on the displacement S and its derivatives is there 
explicitly shown. 



Appendix B 

In this appendix it is demonstrated that the separable second-order solution _E(r)S^ 2 )(q) 
is the most general solution of the Lagrangian fluid equations, once only the second-order 
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terms are retained. Let us suppose that, ab absurdo, the second-order solution is non- 
separable (and longitudinal for simplicity), namely D{r) 2 V$ 2 (q, t) : the factorization of 
the term D 2 does not alter the demonstration. The first order solution is given in Section 
4.1. From the Lagrangian Poisson equation (18), one gets to second-order 

D 2 d 2 ® 2 + 4DDd T §2 + (2D 2 + 2DD - a(r)D 2 ) $ 2 = -a(T)D 2 V 2 , (61) 

where the tilde "~" indicates the Fourier transformed quantities (see next appendix). The 
function P 2 , whose explicit expression is superfluous to give here (the interested reader 
is addressed to the next appendix), is defined by the relation V 2 = — p~ 2 /U 2 (S( 1 )): the 
important point is to note that Vi depends only on the spatial variable: V 2 = ^(p), where 
p indicates the comoving Lagrangian wavevector, and p = |p|. Now it is easy to verify that, 
if the function B{r) satisfying the differential equation K 2 B + 4KB + (2 + aK 2 )B = 1, 
where K = D/D, is introduced, the function $2 m &y be recast in the separable form 

$ 2 ( P ,t) = [2£(t)-1]7Mp). (62) 

Then the function E so defined, E(t) = D(t) 2 [2B(t) — 1], satisfies the differential equation 
E — a(r)E = — a(r)D(r) 2 , which is the first equation in (31). This concludes the demon- 
stration. Similar considerations may be extended to higher-order modes (Buchert, private 
communication; Ehlers & Buchert, 1995, in preparation). 

Appendix C 

In this third appendix, we perform the complete Fourier analysis of the Lagrangian mo- 
tion described in the main text. The final results are useful for practical and numerical 
applications. 

Let us indicate by p the comoving Lagrangian wavevector. The nth-order displacement 
potential (q) may be written as a Fourier integral, 

^ (n) (q) = ^3 / d P V> } (pK p - q , (63) 
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where it is understood that we restrict, in this appendix, to the cases n — 1, 2, 3. Observing 
that ipaW n \p) — S^Hp), h is immediate to obtain from the solution (33) that 

¥ 2) (P) = ~j ^0[{2vf5 D { Vl + p 2 -p)] ^( Pl , P2 )^(Pi)^ (1) (P2) , (64) 
where we have defined the kernel 

k {2 \puP2) = ^ [pIpI ~ (Pi 'Pa) 2 ] = 2 (piP2sin^ 12 ) 2 , (65) 

being e.g. p = |p| and 12 = arcos(pi • p 2 /pi£> 2 ) the angle between the vectors pi and p 2 ; 
the presence of the Dirac-function Sd comes from momentum conservation in Fourier space. 
An alternative simpler expression of ip^(p) may be obtained performing one integration: 
the result may be written as follows, 

¥ 2 \P) = ~j ^3 « (2) (p,p') ^ (1) (P')^ (1) (P-P') • (66) 

We stress that to obtain the expression of in terms of the first-order Zel'dovich potential 
it is completely unnecessary, for practical uses, to specify the form of the divergence- 
free vector R( 2 ), for any realistic initial conditions. The kernel «;^( P i,p 2 ) describes in the 
Lagrangian Fourier space the effects of the non-linear dynamics. 

Similar considerations and calculations may be easily extended to the third-order solu- 
tions. We give here the explicit expressions of the Fourier components of the Lagrangian 
potentials ip^ and originating the longitudinal motion, and of the vector potential 
A, describing the transversal motion. 

To show how the calculations progress, it is more simple to start with the derivation of 
the Fourier component ip^ ' . From the expression (44), 
-(3)/ x_ Pa f rfpi<ip2 



$ S) (P) = ^^[(^MPi + P2 - P)] (flPla - Pi • P2P2«) ^ (1) (Pl)^ (2) (P2) . 

P 71 (67) 

Inserting the solution (66), we eventually obtain, in terms of the Zel'dovich potential, the 
expression 

$ s) (p) = / d j0 4 3) (pi,p 2 ;p) ^ (1) (pi)^ (1) (p2-pi)^ {1) (p-P2) , (68) 

28 



(3) 

where the kernel n b has been introduced: 

4 3) (pi,P2;p) = ^ (2) (pi>P2) 



l 



l - 



P P2 

. VV2 , 



In a similar fashion, the Fourier component ip^ may be derived: 



^ 3) (P) = ~j 4 3) (pi,p 2 ;p) ^ { %i)^ (1 W 



Pi - P2) , 



where, again, we have defined a kernel according to the definition 



4 3) (pi>P2;p) = 7 



6 



€p V a (P/3 - Pl/3 - P2p) Pl v P2a 



(69) 



(70) 



(71) 



e«7<5 Pa PI7 P26 

The algebraic relation S^ C = |e Q7< 5 S^j S^/J has to be used. Finally, the vortical 
component T a along the ^-direction is given by the integral 

Ta(p) =</ 4 3) (Pi,P 2 ;p) ^ {1) (Pi)^ {1) (P2-Pi)^ (1) (p-P 2 ) , (72) 

where the kernel which indeed depends upon the direction, is explicitly given by 

1 



,(3) 



(pi,p 2 ;p) 



,2 " Pl , P2 ) l-^f (2p,-p). 



(73) 



Specifically, the vortical vector T(p) has only two independent components, in that it 
satisfies the condition of transversality, p • T = 0. Therefore it follows that the vector 
potential A is fully specified by only two of its three components: without loss of generality 
we can decide that A z = 0. The significant Fourier components A x and A y are thus given 
by the relations: 

Pz A y (p)=tf x (p), (74) 

and 

p z A x (p)=if y (p). (75) 

These last two equations complete the analysis of the Lagrangian motion up to the third- 
order perturbative approximation: specifically, the third-order dynamics is fully described 
in terms of the tetra-potential \I/ = ipl 3 \ Ai, A2). The complications due to the non- 
linear evolution are summarized in the five kernels k^ 2 \ and the two last ones 
depending on the chosen direction. 
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Appendix D 



In this last appendix, and as an example of how our formalism works, we want to show 
that the general irrotationality condition in Lagrangian space 



£af3y Sf3 a (1 + V • S) 5 ya — S~ /a + S — , 



(76) 



leads to the eq.(46) for the third-order transverse component T of the displacement S. 
From the ansatz (45) one immediately obtains 



Sf* = D{r) Sfl + E{r) Sf> + F c {t) T Pa , 



(77) 



(the F a -mode and i^-mode do not enter in the vortical couplings) and, from eq.(61), we 
get 

= — F c e Q/ 3 7 T 7/3 + 
D^(l) 5, a + £/n(2) 5 7CT - DS$ - ES$ + D 2 S$ C ] , (78) 



DSfJ + ESfl 



i.e. 



e Q/37 {DD'S^ S$ c + DES$> [^(2) 5 7<T - S$] + DES$ [^(1) 5 1(T - Sg>] - F c T 7/3 } = , 

(79) 

where for brevity we wrote /ii(n) = /ii(S^). Recalling that e a /3 y S^ = 0, e a p^S^ = and 
e a/3-yS9'}S^ 1 J c = 0, we finally obtain 



-a/37 



(DE - DE)S$S<$ - F c Tj = , 



JO) 



which coincides with the relation (46) in the main text. 
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